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In central Pb-Pb collisions at LHC the multiplicity of particles produced is so high that it should 
become feasible to examine the nature of transition from the deconfined quark-gluon state to the 
confined hadron state by methods that rely on the availability of high multiplicity events. We 
consider four classes of the transition process ranging from critical behavior to totally random 
behavior, depending on whether or not there is clustering of quarks and on whether or not there is 
contraction of dense clusters due to confinement. Fluctuations of bin multiplicities in each event are 
quantified, and then the event-by-event fluctuations of spatial patterns are analyzed. A sequence 
of measures are proposed and are shown to be effective in capturing the essence of the differences 
among the classes of simulated events. It is demonstrated that a specific index has a low value for 
critical transition but a larger value if the hadronization process is random. 

PACS numbers: 24.60.Ky, 25.75.-q, 25.75.Nq 



I. INTRODUCTION 

Phase transition has always been a subject of great 
interest in many fields. The possibility of observing ev- 
idences for the critical point in heavy-ion colhsions has 
invigorated extensive experimental programs at various 
laboratories [ll, The physics of QCD critical point 
concerns a dense system of strongly-interacting matter 
that is in thermal and chemical equihbrium and is at the 
end of tlie phase boundary between quark and hadron 
phases Q. The search for signals of that boundary in- 
volves experiments at energies where high baryon density 
can be produced. At the Large Hadron Cohider (LHC) 
where the collision energy is much higher, the physics of 
hadronization involves issues that are different from that 
probed at lower energies where high chemical potential is 
expected. Since a hot and dense plasma is produced at 
LHC, the deconfined state persistes for a long time before 
the system is dilute enough to undergo transition to the 
hadron phase. That transition may or may not be rec- 
ognizable as a critical phenomenon, since hadronization 
takes place on the surface over a period of time while the 
system expands. The accumulation of hadrons emitted 
over that period can smear out any signal of interest even 
in the best circumstance for critical transition, which is 
not assured on theoretical grounds. Nevertheless, or per- 
haps particularly because of the difficulty in detecting 
revealing signals, it is of interest to investigate whether 
appropriate measures exist. It is the aim of this work to 
find the most effective way to extract dynamical signals 
indicative of collective behavior of a system transition- 
ing from quarks to hadrons. Our search does not rely 
on the validity of any theoretical view on the possibility 
of critical behavior in heavy-ion collisions at very high 
energy. 

Despite the fact that a large amount of data on Pb-Pb 
collisions at LHC has been produced in the past year Q , 



so far no spectacular departure from earlier expectations 
@ by extrapolation from the Relativistic Heavy-Ion Col- 
lider (RHIC) has been reported. It is not certain whether 
it is because there is no unexpected new physics, or be- 
cause unanticipated physics is not revealed in the con- 
ventional observables. We venture here to ask whether 
there may be an area of investigation of the LHC data 
that has thus far been neglected but may prove to be 
fruitful when suitable measures are used. 

An obviously outstanding feature of the data obtained 
at LHC is that the total multiplicity A'ch of charged par- 
ticles produced is unprecedentedly high, around 6 x 10"^ 
in central collisions at 2.76 TeV Observables that rely 
on large A'ch can thus be exploited in ways that could not 
have been possible at lower energies, thereby forming a 
frontier that has not been explored. Another aspect of 
central Pb-Pb collisions is that there is surely a decon- 
fined system of quarks and gluons, which must undergo 
some kind of transition to the confined state of hadrons. 
Theoretically, the use of Cooper-Frye scheme to cal- 
culate the hadronic properties avoids the issue about the 
nature of the deconfinement-to-confinement transition. 
In our investigation we leave open the question of what 
the nature of that transition is. We propose four sce- 
narios that range from critical behavior on one end to 
non-critical on the other. The focus is on the discovery 
of measures that can possibly distinguish those cases. To 
test whether the proposed measures are feasible experi- 
mentally, we simulate the events in each case and examine 
the effectiveness of the measures. 

The theme of our investigation is about the fluctua- 
tions of spatial patterns during the quark-hadron transi- 
tion. The momentum of each particle is usually expressed 
in terms of {pr, where pr is the transverse momen- 
tum, 1] the pseudo-rapidity, and (p the azimuthal angle. 
For any interval of pr one can study the (r?, 4>) distribu- 
tion in any given event, called the lego plot when shown 
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in some finite-size binning. For convenience, we refer to 
that distribution as a spatial pattern, which changes from 
event to event. The basic question is whether the fluctua- 
tions of those patterns contain any information about the 
nature of the quark-hadron transition. To study those 
patterns one needs good resolution from the experimen- 
tal side and efficient description from the theoretical side. 
Furthermore, one should not integrate over allpT because 
the superposition of different patterns at different ApT 
intervals can smear out all recognizable features. Thus to 
be able to have high resolution in all (px, rj, 0) variables 
as well as to have enough particles in small bins in each 
variable requires high event multiplicities. That is where 
the LHC data become extremely useful. 

There are recent reports from ALICE on the event-by- 
event fluctuations of global observables in Pb-Pb colli- 
sions [1, Q . Non-statistical fluctuations of mean pT are 
found to be lower than a scaling behavior in event mul- 
tiplicity as the collisions become more central. No ex- 
planation of the phenomenon is given in terms of known 
models. Our study here is on local multiplicity fluctua- 
tions in precisely those central collisions. 

There are two kinds of fluctuations that we shall con- 
sider. One is the fluctuation of bin-multiplicity from bin 
to bin in an event, and the other is the event-to-event 
fluctuation of the spatial patterns. On the former it can 
be just random fluctuations, but it can also be in the form 
of clusters of all sizes, as in second-order phase- transition 
[l^ . Both types of patterns will be generated to initiate 
quark conflgurations in the {rj, (j)) space before and dur- 
ing time evolution throughout the period of hadroniza- 
tion for two distinct classes of models. We shall construct 
an algorithm for simulating the global effect of color con- 
flnement through contractions of dense clusters between 
time steps in the hadronization period and the opposite 
effect of thermal agitation by randomization right after 
each contraction. Pion formation according to chosen 
criteria can occur throughout the process, resulting in an 
event distribution of pions in 0). That distribu- 

tion is then analyzed by factorial moments to fllter out 
statistical fluctuations. Repeating the simulation over 
many times generates event-by-event fluctuations which 
turn out to be crucial to the flnding of revealing signa- 
tures. Large deviations from the average event structure 
are possible but very rare; however, they can make sig- 
niflcant influence on the measure that we shall propose. 



II. LOCAL MULTIPLICITY FLUCTUATIONS 
WITH CRITICAL CLUSTERING 

In search for ways to simulate conflgurations that have 
a wide range of characteristics without theoretical preju- 
dices, ranging from critical to non-critical cases, we focus 
in this section on finding a simple procedure to generate 
configurations that correspond to critical behavior. We 
first review a number of related investigations on the local 
behaviors of multiplicity fiuctuations of systems under- 



going a second-order phase transition describable by the 
Ginzburg-Landau (GL) theory. They are then to be con- 
nected to a distribution of cluster production that can 
readily be used to simulate initial configurations before 
hadronization begins. Only one parameter is needed to 
describe the clustering. Later, for the non-critical case 
it is only necessary to change that parameter so as to 
obtain random distributions. 

It was suggested that the normalized factorial moments 
Fq can be used as a quantitative measure of local fluctu- 
ations jlH. Fq is defined by 



where n is the multiplicity in a bin of size S'^ in a d- 
dimensional phase space, and the averages are performed 
over many events. A power-law behavior 

Fq{S)o,6-^- (2) 

over a range of small 6 is referred to as intermittency, and 
has been observed in many systems of collisions [l^l ■ The 
virtue of Fq is that it filters out statistical fluctuations, 
as can be seen as follows. If the multiplicity distribu- 
tion Pn can be written as a convolution of the Poisson 
distribution P^(n) and a dynamical distribution D{ra), 
i.e., 

1"°° to" 

Pn - / dm—e-'^Dim), (3) 
Jo n\ 

then the numerator of Fq is 

V r^^^Pn = / dmm'^Dim), (4) 

which is a simple moment of D{m). Thus, if the dynamics 
is trivial, i.e., D{m) = 5{m — n), then Fq = 1 for all q. 
Any deviation from 1 is then a measure of non-trivial 
dynamical fluctuations, and a power-law behavior in Eq. 
^ would suggest a dynamics that is not characterized 
by a particular scale . The footprint of a PT that has 
fluctuations of all scales may then be possibly observed 
in the measurement of intermittency ^13j . 

To obtain a theoretical quantiflcation of PT in terms of 
Fq , a study of second-order PT in the Ginzburg-Landau 
theory was carried out in Ref. [31, in which the order 
parameter is identified with the multiplicity density. It 
is found that to a high degree of accuracy Fq satisfies the 
power-law behavior 

Fq^F^\ (5) 

where 

Pq = iq- l)^ u = 1.304, (6) 

essentially independent of the details of the GL param- 
eters. Such a behavior was experimentally verified by 
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the study of photon number fluctuations of a single-mode 
laser at the threshold of lasing ^5|| , confirming that it is a 
PT problem describable by GL theory ■ On the theo- 
retical side, it has also been found that using the 2D Ising 
model to simulate quark-hadron PT the resulting scaling 
behavior of Fq is in agreement with Eqs. ([5|) and ([6]) (l7l |. 
It does not mean, however, that an analysis for Fq in the 
current data from heavy-ion collisions can verify or falsify 
the connection between hadronization and second-order 
PT because of the complications that are present in such 
systems but absent in the optical system. The following 
sections are aimed at addressing such complications. We 
note that in the discussion above concerning Eqs. (H]) to 
Q no mention is made of the specifics about multiplici- 
ties, whether hadrons or quarks. In the remainder of this 
section we shall use particles as a generic term that can 
refer to either hadrons or quarks. The result to be ob- 
tained will be used in the following sections to generate 
configurations of quarks (and/or antiquarks) just before 
hadronization. 

From a different perspective it is of interest to study 
the problem of clustering, since critical phenomenon is 
known to be characterized by clusters of all sizes. It is a 
subject relevant to our investigation here because we shall 
analyze the 2D spatial distribution of particles in the rj-cj) 
space. Cluster formation in the context of hadronization 
has been considered before but from the point of 
view of self-organized criticality instead of second-order 
PT. Since the latter exhibits the behavior of Eqs. (O and 
([6]), it is not difficult to see how such a behavior can be 
achieved by an adjustment of the scaling property of the 
cluster formation. 

As a generic problem on clusters of particles, let C be 
the number of particles in a cluster, and P{C) be the 
probability distribution in C, which we assume to have 
the scaling form 

P{C) oc C-^, 7 > 0. (7) 

We let the center of a cluster be distributed randomly 
in a 2D space which we take to be 1 unit of length on 
each side. We further let a particle in a cluster C to be 
distributed randomly around its center with a Gaussian 
width of 

a = 0.1 (8) 

so that there can be high density of particles, though 
with decreasing probability at high C. Such spikes of 
multiplicity are what can give rise to intermittency. 

To investigate the multiplicity fluctuations in small bin 
sizes, we divide the unit square into AI x M bins, M vary- 
ing from 8 to 70. We let iVo particles be distributed in the 
unit square in accordance to Eq. ([7]) and (jS]), and allow 
A^o to be approximately 100, since the sum of all clusters 
may not be exactly the same for every event. With n 
being the number of particles in a bin, we perform the 
calculation for Fq{S) according to Eq. by averaging 
over all bins first, and then averaging over all (5 x 10^) 



events for any fixed 5 = 1/M. Only bin multiplicities 
with n > q are counted in Fq(M). The result for 7 = 2 is 
shown in Fig. 1(a) that exhibits non-trivial dependence 
on M; it is not strictly linear in the log-log plot. How- 
ever, when Fq is plotted against F2 for q — 3,4,5, very 
good linearity is found, as shown in Fig. 1(b) replicat- 
ing the behavior found for GL [l^ . Similar power-law 
behavior is obtained for other values of 7. We show in 
Fig. 1(c) the dependence of the exponent (3q, defined in 
Eq. ([5]), on g — 1, and find in accordance to Eq. © that 
I' = 1.136,1.171, and 1.315, for 7 = 1.5,2.0,2.5, respec- 
tively. 
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FIG. 1; Intermittency analysis of clustering model for 7 = 2 
in Eq. ((T]): (a) scaling in M, (b) _F-scaling in Fq vs F2, and 
(c) power-law behavior of Pq in Eq. ((6]). 



What we have done above is to find a quick way to 
relate a particle distribution in a 2D space to a GL-type 
-Fq-scaling behavior through the use of a cluster distri- 
bution. Since a range of values of v has been obtained, 
the relationship between 7 and v has no specific dynam- 
ical significance. It merely demonstrates that large local 
multiplicity fluctuations can generate large Fq with scal- 
ing properties. In the case of 7 = 2.5 we see in Fig. 1(c) 
that the straight line connecting the three points, when 
extrapolated, does not go through = 1, thus showing 
slight deviation from Eq. ([5]). The special value v = 1.3 
for second-order PT can approximately be achieved by 
using a value of 7 between 2.0 and 2.5. Although no 
dynamics has been introduced to establish that connec- 
tion, our objective here is accomplished by having found 
a simple procedure to simulate configurations that may 
be relevant to a quark and antiquark system that is at the 
edge of PT to hadrons. This is only the starting point of 
a much more complicated problem in heavy-ion collisions 
that we shall describe next. 
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III. SIMULATION OF EVENTS WITH SPATIAL 
FLUCTUATIONS IN HEAVY-ION COLLISIONS 

If there is no dynamical structure in the spatial pat- 
tern of soft hadrons in the lego plot (of 7] and (p), then 
there should be nothing interesting in whatever measure 
that is used to analyze the data. The reverse is, however, 
not necessarily true. Unless a measure is sensitive to the 
consequence of certain dynamics, one may not observe 
what is interesting. To provide motivation for an experi- 
mental effort to search for an unconventional signal, it is 
necessary to demonstrate the worthiness of such an effort 
under the best of circumstances. The optimal scenario 
is that the quark-gluon system undergoes a second-order 
PT in forming hadrons so that there are large multiplicity 
fluctuations. It is not known whether such a phenomenon 
occurs at LHC. Our aim here is to simulate events that 
belong to different classes of dynamical characteristics, 
ranging from robust criticality to mundane randomness. 
For each class of hadronization the measure to be used 
in the following section should exhibit distinguishing fea- 
tures in the hope that analysts of the real data would 
have the incentive to pursue the more difficult task of 
extracting worthwhile information from what is actually 
observed. 

The physics we aim to simulate is for the duration be- 
tween the end of the quark phase and the beginning of 
the hadron phase. It is assumed that the density is so 
low that the confining forces among the quarks and the 
antiquarks begin to redistribute them spatially. Gluons 
are assumed to have been converted to q and q that are 
the basic units prior to hadronization. We restrict our 
attention to the central rapidity region with |?7| < 1 and 
< < 27r for central collisions. The surface of the 
plasma cylinder with those values of rj and (j) is mapped 
to the unit square S. Since the fluctuation properties 
to be discussed below do not depend on the precise area 
that S corresponds to, a portion of the (?], 0) region in 
the actual data should suffice to serve as a workable ba- 
sis for analysis. The cylinder hadronizes at the surface 
only, layer by layer. Thus we have to consider many time 
steps ti,i = 1,2,---. Between adjacent time steps the 
q and q adjust their positions in S in accordance to a 
procedure that we impose to simulate conflnement and 
pion emission. The quarks have thermal px distribution 
whose inverse slope Ti decreases incrementally with ti. 
At a new time step a new set of q and q are introduced 
to represent the movement of the next layer of quarks to 
the surface. The conflnement procedure starts over again 
until most of the quarks (not necessarily all) are pionized 
before another time step is taken. This is a general out- 
line of the algorithm to convert q and q to pious, whose 
coordinates in (px-, r], (j)) are registered for event-by-event 
analysis later. 

The main part of our physics input is the deconfine- 
ment to conflnement transition between time steps. The 
principal characteristic of a critical phenomenon is the 
tension between the ordered and the disordered motions 



of a system. In the Ising model of a magnetic system the 
near-neighbor interaction between spins tends to align 
all spins in the same direction, but the thermal motion 
tends to randomize them. If the quark-hadron PT is in 
the same class of criticality, then the tension is between 
conflnement that draws the qq into an ordered pair and 
the disordered thermal agitation that tends to keep the 
quark system in the deconflned state. On the surface 
of the plasma cylinder the q and q may start to cluster 
as the density and temperature get close to the transi- 
tion point, thereupon a local region of high density con- 
tracts under the conflnement force to improve the likeli- 
hood of a qq pair to fall within a conflnement distance. 
In principle, we should keep track of the colors of the 
quarks, but that would complicate the procedure much 
more. The simple algorithm we adopt mimics the general 
idea of color confinement and generates large multiplic- 
ity fluctuations. One of the chief objectives of this work 
is to see whether many such local fluctuations at differ- 
ent time steps can survive the superposition at the end 
of the hadronization process, and still be detected by an 
effective measure. The feasibility issue is addressed at 
the cost of precision in QCD, which is difficult in the soft 
regime over an extended period. 

Let us describe first the case of critical hadronization, 
followed later by other cases that are less critical. 



A. Critical 

1. Initial Configuration 

We start by seeding the unit square S with 500 pairs 
of qq. They are clustered according to the probability 
distribution P{C), given in Eq. d?]). That is, C pairs of 
qij are placed in a cluster centered at a random point 
in 5*; they are Gaussian distributed around the center 
with a width specified by Eq. ([5]). As a result of the 
study in Sec. II that relates the exponents 7 and v, we 
choose 7 = 2, corresponding to v not quite up to 1.3. 
As the cluster multiplicity C is summed over all clusters, 
we stop the seeding process when the total just exceeds 
500. Inside each cluster the q and q are not correlated; 
they are independently distributed in (77, 0) as well as in 
Pt, for which the thermal distribution is the exponential, 
exp(-pT/T), with T set at 0.4 GeV. 

The above procedure is for setting up an initial con- 
figuration at to = 0, counting from just before the first 
pions are emitted, but long after the collision time. The 
value of T is not set at Tc which is lower, but at a value 
above the observed average {pt), since the bulk of the pi- 
ons will be produced at later time when T will be lower. 
The clustering is put in to generate a spatial configu- 
ration that is most likely to represent a system moving 
toward a critical transition. 
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2. Pionization 

In the above configuration if a g and q are within a dis- 
tance d from each other that is less than do = 0.03, then 
we regard the pair as effectively a pion and let it be taken 
away from the configuration, but registering the pion po- 
sition in a separate (77, 0) space at the midpoint between 
the pair. We assign a value of pr to that pion that is 
equal to the sum of the values of the q and q. This is 
based on the recombination model, where the recombi- 
nation function has a momentum-conserving (5-function: 
5{pt — Pit — P2t) [T9l - [2lj . The thermal distributions of 
quarks and pions have the same T . We ignore the color 
and flavor of the quarks without losing the essence of lo- 
cal fluctuation, since if the q and q had color and flavor 
labels, it would take longer in the iteration process for 
the qq pair in a cluster to pionize without changing the 
(77, 0) coordinates appreciably. 



3. Contraction 

The probability that a qq pair is within dp apart is 
small, so the majority of the quarks remain in S. They 
are under the influence of the color forces to move toward 
confinement. To describe that movement in a collective 
way rather than in terms of pair- wise interaction that is 
unrealistic, we adopt the contraction procedure as fol- 
lows. 

Let S be divided into 5x5 bins. Calculate the num- 
ber of q and q in each bin. Separate the bins into two 
types: dense bins have more q and q than the average 
bin multiplicity, and the dilute bins have less. The dif- 
ference between q and q is ignored here. If adjacent dense 
bins share a common side, they are grouped together as 
members of a cluster of dense bins. Let D refer to such a 
cluster of dense bins, which may have an irregular shape, 
but are connected. Let No denote the number of bins in 
D. Define ro to be the coordinates in S that is the center 
of mass of D. Now, we do a contraction of D. That is, 
we redistribute all the q and g in Z), centered at but 
with a Gaussian width 

(TD = (7iN'J\ (9) 

where tJi is a parameter that characterizes the degree of 
contraction. We use cri =0.1 here with other possibilities 
to be discussed later. Because of the Gaussian distribu- 
tion, the q and q that are spread out originally in Nu bins 
are drawn closer together to be located mostly within the 
Gaussian peak, resulting in a contraction. Since there are 
5x5 bins, the original bin size is 0.2 x 0.2, to which an 
should be compared after contraction. This is how we 
model the effect of the ordered motion due to confine- 
ment. We repeat step 2 to allow qq pairs to pionize in 
this new configuration. 



4- Randomization 

The disordered motion that counter-balances the or- 
dered motion is the thermal randomization. We model 
that by requiring all the q and q in the dilute bins to be 
redistributed randomly throughout S, resulting in a new 
configuration. We then repeat step 3 to have another 
round of contraction and pionization. Each time in the 
iteration process more qq pairs are converted to pions, as 
we alternate contraction and randomization until around 
95% of the qq system is depleted. We regard that as the 
end of one time step in which the quarks on the cylin- 
der surface are hadronized. We then proceed to the next 
time step when a new layer of quarks move up to the 
surface. 

We note that before the next time step is taken, pi- 
onization becomes increasingly difficult when fewer and 
fewer q and q remain to find their partners to coalesce. 
Decreasing ai to 0.05 speeds up the process near the 
end, but does not change the pion distribution in (77, 0) 
very much. Physically, it is not necessary that all qq in a 
layer hadronize before the next layer of quarks moves up. 
Leaving roughly 5% to be mixed with the next layer of qq 
seems reasonable in our attempt to model a continuous 
process of hadronization by discrete steps. 

5. Subsequent Time Steps 

To advance to the next time step we introduce 100 new 
pairs of qq according to the distribution P{C) and add 
them to the existing q and q that remain from the pre- 
vious step. We then follow the same procedure as above 
to contract, pionize, and randomize repeatedly until 5% 
of q and q remains. At the ith time step all q and q have 
the pt distribution with an inverse slope 

T, = r,_i - 0.02 GeV. (10) 

That is the Ti that the pions emitted at U will also have. 
This is carried out 10 times so that in total we intro- 
duce 1500 qq paris, most of which turn into pions. To 
have approximately 1400 pions produced in jr/j < 1 in 
central collision at LHC is not unrealistic. Since we are 
not interested in global multiplicity fluctuations in this 
study, we stop after ten steps. The pions collected in 
{pt, Tj, (j>) for the event is then stored for later analysis of 
local fluctuations. 



B. Quasi-critical 

Consider now the case where clustering does not occur 
in the initial conflguration, nor when new layers of q and 
q move to the surface. Thus the initial 500 pairs of qq 
are seeded randomly in S, and so are the subsequent 100 
pairs at each time step. There is then no correspondence 
with the second-order PT discussed in Sec. II. Techni- 
cally, this corresponds to the exponent 7 in Eq. ^ being 
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very large, because the probability for having a large clus- 
ter C is then very small. Specifically, we choose 7 = 5 for 
no clustering. One may regard this case as being in corre- 
spondence to a cross-over in the phase diagram where no 
distinct boundary between the quark and hadron phases 
can be identified. However, we require that the confine- 
ment forces to be still at work to turn quarks to hadrons. 
So we use the contraction-randomization procedure de- 
scribed in subsection A to carry out pionization. Such 
a procedure simulates the tension between confinement 
and deconfinement, but cannot be regarded as what is 
needed for a critical transition. The coordinates of each 
pion in [pT, r], <p) are recorded as before for later analysis. 



at the position midpoint between the qq pair. Then we 
go to another quark and repeat the process. When no 
more pairs can be found within that distance, we pro- 
ceed to the next time step whatever numbers of q and q 
are left. With more qq pairs supplied from the next layer, 
the probability of pionization is increased. This process 
will not convert all qq pairs to pions since the quarks can 
be far apart from antiquarks without contraction. To 
have quarks left over at the end of ten time steps does 
not matter as far as the local multiplicity fluctuation is 
concerned. We shall see that there is still interesting 
structure that can be extracted from the accumulated 
events simulated that way. 



C. Pseudo-critical 

Suppose now that we seed the configurations with clus- 
tering, but do not impose contraction between time steps. 
Thus the q and q configurations, as each layer reaches the 
surface, are close to the critical condition, when we set 
7 = 2 as we have done above. However, the hadronization 
process is carried out by letting qq pairs form pions when- 
ever a pair gets close together. Without contraction we 
cannot require the distance between pairs to be less than 
do as in Subsec. A. 2, since the probability for that to oc- 
cur is low. Without new dynamics during hadronization 
we simply let each quark to search for the nearest anti- 
quark within a distance d = 0.1 and let a pion be formed 



D. Non-critical 

To the other extreme situation away from the above 
three cases we consider the non-critical case of no orga- 
nized dynamics at all, i.e., random configurations and no 
contraction. The result should not have any content of 
interest. We have nevertheless carried out the simulation 
and will show the result that is non-trivial and therefore 
instructive. 

A summary of the four cases can be expressed as a 
matrix shown below. 





clustering 


no clustering 


contraction 


critical 


quasi— critical 


no contraction 


pseudo— critical 


non— critical 



(11) 



IV. MOMENTS FOR EVENT-BY-EVENT 
FLUCTUATIONS 

We now consider the method of analysis of the many 
events either as measured at LHC or as generated in 
the models described in the preceding section. Because 
of the high multiplicity of particles produced in cen- 
tral collisions, we make cuts in pT, such as in an in- 
ternal ApT around pt = 1 GeV/c. We have considered 
Apr = 0.04,0.07 and 0.1 GeV/c. In such small intervals 
the particle multiplicities are significantly reduced and 
spatial patterns in (jj, <p) begin to appear because of the 
possibility of empty bins. 

To find an effective measure of the fluctuations we need 
it to be sensitive to both the spatial variations from bin 
to bin and the event-by-event fluctuations. We refer to 
the former as horizontal and the latter as vertical. Differ- 
ent moments are to be taken for horizontal and vertical 



averages so as to allow spatial and event- wise fluctuations 
to manifest separately. 

We divide the unit square into bins with M being 
not more than 70, depending on the multiplicity in the 
ApT interval, just so that the important part of the M 
dependence is captured. For the spatial fluctuations we 
use the horizontal factorial moments for event e, defined 
as 

F,%M) = /|(M)/[/f (M)]", (12) 

where 

f^{M) = {n{n-l)--.{n-q + l))^. (13) 

The average ( )/i is performed over all bins, n be- 
ing the multiplicity in a bin. Only n > q is counted in 

/|(M). Clearly, if M is very large, fq(M) may be zero, 
for q>2, although /f (M) is never zero. However, there 
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may be an event where (M) may not vanish at a large 
M; then it would imply sharp spikes of multiplicity in 
some bins. If the fluctuations among the bins are Poisso- 
nian, then we have F^{M) = 1 for any M, as is the case 
with vertical fluctuations discussed in Sec. II. Interesting 
spatial fluctuations are, however, not Poissonian. 

For each event we can calculate F^{M). We note that 
F^{M) is a simple characterization of spatial pattern, 
but is not the only one possible. Any alternative descrip- 
tion can also be used in the study below on event-by- 
event fluctuations. One may therefore regard Fq{M) as 
a generic symbol. 

If {Fq{M))^ denotes the vertical average of Fq{M) over 
all events, then the fluctuation of i^|(M) from {Fq{M))^ 
is what we want to quantify. To that end we consider the 
pth-order moments of 



deconflnement-to-conflnement transition, it will become 
clear in the next section, where model calculations can 
generate concrete local multiplicity fluctuations, that it 
is more comprehensive to consider the range 1 < p < 2 in 
the study of Cp^q{M). If in that range V'g(p) has a linear 
dependence on p, it is better to define the slope in that 
wider range as 



pi, = dil}q{p)/dp. 



(19) 



%{M)^FI{M)I {Fq{M))^. 



(14) 



It is this quantity ^q, referred to as erraticity indices, 
that will become an effective measure of the criticality 
classes that is independent of p and M . A large value 
of /ig means that there are influential contributions to 
Cp,q{M) from large Fq{M) weighted more heavily at 
large p, which in turn implies that very erratic fluctua- 
tions of the spatial patterns are involved to render Fq{M) 
non- vanishing at large q and M . 



I.e. 



Cp,,(M) = (cI>^(M))^, 



(15) 



which is a double moment introduced earlier for the study 
of chaotic behavior of particle production in branching 
processes p^ . Whereas q must be an integer, p need not 
be. In fact, the derivative of Cp^q{M) at p = 1, i.e.. 



I]g(M) = -i^Cp,q{M) 



11 V 



(16) 



has been related to an entropy defined in the event space 
[2^ , but by itself it is not very useful because it depends 
on M . Simplification can occur if Cp_q(M) has a power- 
law behavior in M 



Then an (entropy) index can be defined as 



„(i) = ^^'?(^'^) = ±^ („) 
^« dlnM dp^'^^P' 



(17) 



(18) 



p=i 



,(1) 



which is independent of M. It was found that /ig ' can 
characterize the fluctuations of spatial patterns so well 
that they are as useful as the Lyapunov exponents in 
providing a quantitative measure of classical chaos [23| . 

The power-law behavior of Cr, .n(M) in Eq. ([T7| has 
been referred to as erraticity [2J, . It was proposed 
as the next logical step to take beyond the intermittency 
analysis. Attempts have been made to find experimen- 
tally the erratic fiuctuations of Fq from event to event in 
multiparticle production fl^ . In meson-proton collisions 
at 250 GeV/c beam momentum, it was found that the 
erracticity measures are dominated by statistical fluctu- 
ations at such low energy [2^. In nucleus-nucleus colli- 
sions interesting signals have been found both in mod- 
els [131 and in emulsion experiments [l^ even at low 
energies. At high energy where we can examine the 












FIG. 2: Examples of bin multiplicity fluctuations in (rj, (j)) for 
the four cases arranged in the matrix form of i.e., (a) 

critical, (b) quasi-critical, (c) pseudo-critical, (d) non-critical. 



V. RESULTS OF MODEL CALCULATIONS 

We have simulated in the order of 10® events for the 
four models ranging from critical to non-critical cases 
described in Sec. III. To give visual images of their qual- 
itative differences, we select one event from each case to 
illustrate their behaviors in the lego plots. For an event 
to contribute to at g > 2, it cannot be typical if the 
average bin multiplicity (n) is small. We show untypi- 
cal events that have large n in some bins and thus can 
contribute to non-trivial F3 at M = 30 and ApT = 0.1 
around pt — ^ GeV/c, for which (n) is approximately 
0.1. In Fig. 2 the four lego plots are arranged as a ma- 
trix in the classiflcation given in ((TT|) . Each of the plots 
corresponds to an event in a given class that contribute to 
the largest value of F3(M) at M = 30; thus there should 
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be at least one bin that has a bin multiplicity n > 3. 
It is evident that the distinctive difference among them 
is that the pseudo-critical and non-critical cases in the 
lower panels (c) and (d) have particles distributed more 
widely throughout the base square than in the other two 
cases above in (a) and (b), where a localized cluster sur- 
rounds a high peak. The intent of this figure is only to 
give a sense of the qualitative difference in the nature of 
spatial patterns when a stringent demand is placed on 
the bin multiplicity to depart from the average. 
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FIG. 3: Intermittency analysis for the critical case. Panels 
are the same as in Fig. 1. 
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FIG. 4: Intermittency analysis for the quasi-critical case. 
Panels are the same as in Fig. 3. 

To see the scaling behavior, we examine the vertical 
averages {Fq{M))^ vs M in log- log plots for g = 2, • • • ,5 
and for the cut Apj- = 0.1 GeV/c. Hereafter the sub- 
script V will be omitted for brevity. For the critical case 
the results are shown in Fig. 3 where in (a) there is an 
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FIG. 5: Intermittency analysis for the pseudo-critical case. 
Panels are the same as in Fig. 3. 
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FIG. 6: Intermittency analysis for the non-critical case. Pan- 
els are the same as in Fig. 3. 



increase with M before M gets larger than 20, but they 
all have similar behavior for different q, so when {Fq) is 
plotted against (^2)1 we see in (b) a simple straight-line 
behavior for q = 3,4,5. Using Eqs. ^ and ^ to de- 
scribe the power jiq, we find in (c) that Vc-cit = 1.41. 
Similar properties are found in the quasi-critical case 
with ^'quasi = 1.33, as shown in Fig. 4. In the pseudo- 
critical cases, shown in Fig. 5, we see robust scaling be- 
havior in (a). F-scaling behavior in (b) yields the value 
t'psoudo = 1-26 shown in (c). Finally, in the non-scaling 
case we find the opposite situation where {Fq{M)) de- 
creases with increasing M, as exhibited in Fig. 6 (a). It 
means that bin multiplicities do not get large enough de- 
viations from (n) through random fiuctuations so that, 
when the bin size gets small, {Fq{M)) approaches 1 that 
one expects from Poissonian fiuctuations, as stated just 
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below Eq. dH). In (b) we still see regularity in (Fg) plot- 
ted against {F2), but it is important to recognize that 
the low end corresponds to high M with all (Fg) around 
1, while the high end is for low M, quite contrary to 
the three cases in Figs. 3-5. One can extract the values 
of (3{q) as in (c), but there is no sensible value of to 
be assigned to this case. In short, the non-critical case 
does not yield any interesting result from the study of 
this kind. We have to go beyond simple intermittency 
analysis in order to find a suitable description that can 
render quantitative comparison between the critical and 
non-critical cases. 

Although the results shown in Figs. 3-5 are clear and 
easily quantifiable in terms of i/, a great deal of infor- 
mation is lost by calculating those averages. To exhibit 
the degree of fluctuations of the spatial patterns from 
the average, let us use P{^g) to denote the probability 
distribution of event ^g{M) at fixed M, where ^g{M) is 
defined in Eq. dHl). In Fig. 7 we show P{^g) for q = 2 
and 4, and for clarity only for Af = 8 and 30. The four 
classes of criticality are again in the matrix format of 
PT|) . We see that in all cases the distributions for g = 2 
and M = 8 (solid lines) are peaked at $g = 1. But for 
other values of q and M, the four cases differ in differ- 
ent ways. In the other extreme situation corresponding 
to g = 4 and M = 30 (lines with crosses), we see that 
P($4) is peaked at $4 = in all cases. That is because 
in small bins the average bin multiplicity (n) is much less 
than 4, so the values of $4 for many events are 0. In fact, 
for (c) pseudo-critcial and (d) non-critical, only a small 
fraction of events have large enough bin fluctuations to 
render F^ non-zero, so P{^4) has a 5 function peak at 
<I'4 = 0, whose positions are shifted in Fig. 7 for visibil- 
ity's sake. In the intermediary values oi q — 4, M = 8, 
the dashed lines in all four cases are all very broad, signi- 
fying wide fluctuations. For g = 2, M = 30 (dash-dotted 
lines) the peaks are around $2 = 1, similar to the solid 
lines. 

In the insets of Fig. 7 (c) and (d) we show that for 
q — 4 and M = 30 there are contributions to P($4) 
at extremely large 4>4 (in the order of 10'^). The verti- 
cal axes have the scale factor 10^^. They balance the S 
functions at $4 = so that the average is ($4),„ = 1, by 
definition. This irregular behavior reveals the nature of 
fluctuations from event to event. When the average bin 
multiplicity (n) is about 0.03, it is difficult to find events 
in which there is a bin with n > 4, unless there are dy- 
namical effects (such as confinement contraction) to in- 
troduce large fluctuations. In cases (c) and (d) almost 
all events have n < 4 so $4 = 0, except for some very 
rare events that make non-trivial contribution to non- 
zero $4, whose values are therefore very large because 
F4{M — 30) is exceedingly large (due to the smallness 
of ff) even though $4 is normalized by {F^^M = 30)) 
which is proportional to the rarity of such events. 

The probability distributions P{^g) contain too much 
information that cannot easily be conveyed. We learn 
from Fig. 7 that for g = 2 there are no drastic differ- 



— q=2,M=8 
- -. q=2,M=30 




---q=4,M=8 
^q=4,M=30 



FIG. 7: Probability distributions of ^q{M) for various values 
of q and M . The panels correspond to the criticality cases 
expressed in (3.3) and as shown in Fig. 2. The insets in (c) 
and (d) are for large values of "I>q as explained in the text. 
The vertical axes have a common scale factor of 10~^. 



ences among the four cases when M is increased from 
8 to 30. It means that bin multiplicities in each case 
can fluctuate sufficiently to exceed n = 2 and generate 
a modest width of the peaks in ^($2) around $2 = 1- 
We therefore should not expect a good measure at q = 2 
to distinguish the criticality classes. For g = 4, however, 
we see significant differences between the cases (a,b) with 
contraction and (c,d) without contraction. To quantify 
their differences we consider the moments Cp^g{M) de- 
fined in Eq. (jlSp. For p = 2, g = 4, we show in Fig. 8 the 
M dependence for three intervals of Apr : 0.04, 0.07 and 
0.1 GeV/c. The cases a, b, c, and d in the legend corre- 
spond to the panels in Fig. 7 arranged in the matrix form 
of (fTTj) . Some of the open symbols are displaced slightly 
from M = 8,16 and 30 in order to avoid them from being 
covered by the filled symbols. We see that they satisfy 
power-law behavior very well in all cases, validating the 
meaningfulness of the erraticity exponents ipg{p) defined 
in Eq. (fT7|) . Similar study can be done for p = 1.25, 1.5 
and 1.75, and yield similar scaling behavior, which we do 
not show for brevity. 

The values of 'ip4{p) are shown in Fig. 9. In that fig- 
ure the cases (a) and (b) are depicted collectively by open 
symbols, while the cases (c) and (d) are by filled symbols. 
The values of ip4{p) in the (c,d) group are all larger than 
those in the (a,b) group. Evidently, there are good linear 
dependencies on p between 1.25 and 2 in all cases. The 
straight lines drawn through them are extended only to 
p = 1.1. At p = 1 it is necessary that ipg{l) = because 
Ci,g{AI) = 1 for any g and M. The lines, if extrapolated 
further to p = 1, would all miss the origin (1,0) by little 
bits. There are some reason for that to happen at the 
point p — 1, which we discuss below. For now, we con- 
centrate on the linear portion in Fig. 9 and determine 
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FIG. 8: Power-law behavior of C2,4(M) exhibiting erraticity. 
The labels a, b, c, d correspond to the four panels in Figs. 
2 and 7 in the matrix format of (lllf) . Pairs of symbols are 
approximated by one straight line, so 12 types of symbols are 
represented by 6 lines. ApT are in the units of GeV/c and 
are around the value pr = 1 GeV/c. 




FIG. 9: Linear dependencies of the erraticity exponents ip4{p) 
on p for the (a,b) cases in open symbols and for the (c,d) cases 
in the solid symbols. Their slopes give the values of the indices 
/X4. The units of Apr are in GeV/c. 



the slopes fj.^, as defined in Eq. ([T^. It is clear that the 
(c,d) group has essentially no dependence on ApT- Even 
in the (a,b) group the spread due to different Apx cuts is 
not severe enough to make it unreasonable for us to give 
the three lines in each group an average slope as follows: 



= 1.87 ±0., 



= 4.65 ±0.06. 



(20) 



It is remarkable that we can obtain numerical summary 
of the different cases independent of p and Af and only 
mildly dependent on Apx- The values of these indices 



mean that for the models (a,b) that have contraction due 
to confinement there are spikes in bin multiplicity locally, 
making them easier to have bins with n > 4 (thus lower 
values of /i4 or less erratic) than for the models (c,d) that 
have no contraction. In the latter cases there are not 
enough multiplicity fluctuations to have n > 4 in most 
bins and in most events, so the probability of having a 
non-zero at high M is extremely low, but in those rare 
events the value of F4 is so high (thus more erratic) that 
the pth moment with p > 1.25 give them such a high 
weight as to raise '04 (p) and /i4 significantly above those 
for (a,b) cases. 

Turning now to the complication at the point p = 1, 
we note that if the straight lines in Fig. 9 go through the 
origin (1, 0) exactly, then the slopes are constant through- 
out 1 < p < 2 and we could alternatively focus on the 
point p = 1 exclusively. In that case we could return 
to Eqs. ((TB| and (IT51) and determine ^l^"^ directly from 
I]4(Af). Unfortunately, we have found that I]4(Af) does 
not depend on £nM linearly, so we are unable to make 
use of the Eq. ([T^ to calculate the slope at p — 1. This 
is consistent with the fact that the straight lines in Fig. 
9 do not cross the origin precisely, and that some bend- 
ing of those lines near the point (1, 0) reveals the lack of 
universality of the slopes for all p. Stated differently, the 
power-law behavior of Eq. (I17|) is not true for all p; the 
local multiplicity fluctuations can be so severe and are 
so different among the four cases that the indices /ig^-* at 
p = 1 cannot summarize their differences. 

Returning to the results expressed in Eq. ((20)) . we can 
conclude that what separates (a,b) from (c,d) is the dom- 
inating effect of contraction over clustering. Recall from 
([TT|) that the two columns are for cases where initial and 
reseeded configurations between time steps are with clus- 
tering (a,c), and without clustering (b,d), while the two 
rows are for cases where the configurations between time 
steps undergo contraction (a,b), and no contraction (c,d). 
Since there are 10 substeps within each time step, con- 
tractions rearrange the configurations sufficiently so that 
the erraticity moments Cp^q retain essentially no mem- 
ory of the clustering effects in the input. Physically, it 
means that as the quark-gluon plasma approaches low 
density near the end of its expansion, whether or not 
there is critical clustering, the confinement forces that 
act on the quarks near the surface exert the dominant 
effect on drawing the qq pairs together in order to pio- 
nize, despite the opposing tendency to randomize due to 
thermal activities that persist. The same process is re- 
peated each time a new layer moves to the surface. The 
tension between confinement and deconfinement is what 
leads to large local fluctuations evidenced by the low val- 
ues of the erraticity indices /i^'^. Without that tension 
the bin-multiplicity fluctuations have no dynamical push 
beyond randomness, so it is highly erratic to have some 
rare events to contribute to non-trivial Cp^4, hence higher 

values of fi^''^. 



11 



VI. CONCLUSION 

The purpose of this work is mainly to describe an un- 
conventional method to analyze the LHC data in the 
hope that some experimentalists may find it adventure- 
some. The new frontier opened up by the high multiplic- 
ity events provides a fertile ground for exploration that is 
not feasible at lower energies. That is why the subject is 
not among those predictions that could be extrapolated 
from RHIC. If there is any hint of critical behavior in the 
quark-hadron transition, that would be a new discovery 
at LHC. Even if nothing critical is found, analysis along 
the line suggested here should lead to deeper understand- 
ing of the hadronization process. 

Our models of the four classes of criticality may not 
turn out to be realistic, but they have been useful in 
testing the effectiveness of the erraticity moments and 
indices. The basic issue is, of course, whether the pro- 
posed measures can be applied to the real data to uncover 
interesting physics. There exist other physical processes 
that are totally ignored in this study. Chief among them 
is minijet production, which is understood to be copious 
in Pb-Pb coUisions at 2.76 TeV i293. Usual jet study at 



LHC is for pr very large, e.g. > 50 GeV/c. Our analysis 
here is for px ~ 1 GeV/c, which is minute compared to 
the jet towers, but the result can nevertheless be affected 
by thermal-shower recombination. It is unlikely that the 
low pt enhancement by minijets can lead to multiplic- 
ities greater than 4 in very small bins, but rare events 
with large fluctuations are what the erraticity moments 
Cp^q{M) are sensitive to. Thus at this point we cannot 
rule out contamination of the phase transition effects by 
minijets. Such possibilities perhaps would encourage the 
experimentalists to investigate the subject, either to find 
resolution of ambiguities or to gain new perspective on 
an aspect of physics that is not well understood. 
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